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RADU BALAN 

Abstract. This paper is concerned with the question of reconstructing a vector in a finite- 
dimensional real or complex Hilbert space when only the magnitudes of the coefficients of 
the vector under a redundant linear map are known. We present new invertibility results as 
well an iterative algorithm that finds the least-square solution and is robust in the presence 
of noise. We analyze its numerical performance by comparing it to two versions of the 
Cramer-Rao lower bound. 



1. Introduction 

This paper is concerned with the question of reconstructing a vector x in a finite-dimensional 
real or complex Hilbert space H when only the magnitudes of the coefficients of the vector 
under a redundant linear map are known. 

Specifically our problem is to reconstruct x G up to a global phase factor from the 
magnitudes , 1 < k < m} where {/i, ■ ■ ■ , fm} is a frame (complete system) for H. 

A previous paper |1] described the importance of this problem to signal processing, in 
particular to the analysis of speech. Of particular interest is the case when the coefficients 
are obtained from a Windowed Fourier Transform (also known as Short-Time Fourier Trans- 
form), or an Undecimated Wavelet Transform (in audio and image signal processing). While 
15 presents some necessary and sufficient conditions for reconstruction, the general problem 
of finding fast /efficient algorithms is still open. In [3] we describe one solution in the case of 
STFT coefficients. 

For vectors in real Hilbert spaces, the reconstruction problem is easily shown to be equiv- 
alent to a combinatorial problem. In [5] this problem is further proved to be equivalent to a 
(nonconvex) optimization problem. 

A different approach (which we called the "algebraic approach") was proposed in [2]. While 
it applies to both real and complex cases, noisless and noisy cases, the approach requires 
solving a linear system of size exponentially in space dimension. The algebraic approach 
mentioned earlier generalizes the approach in ^ where reconstruction is performed with 
complexity O(n^) (plus computation of the principal eigenvector for a matrix of size n). 
However this method requires m = 0{n'^) frame vectors. 

Recently the authors of [12] developed a convex optimization algorithm (PhaseLift) and 
proved its ability to perform exact reconstruction in the absence of noise, as well as its 
stablity under noise conditions. In a separate paper, [13], the authors developed further a 
similar algorithm in the case of windowed DFT transforms. 
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In this paper we analyze an iterative algorithm based on regularized least-square criterion. 
The organization of the paper is as follows. In section |2] we define the problem explicitely. 
In section [3] we describe our approach and prove some convergence results. In section |4] we 
establish the Cramer-Rao lower bound for benchmarking its performance which is analyzed 
in section [5j Section |6] contains conclusions and is followed by references. 

2. Background 

Let us denote by H the n-dimensional Hilbert space M" (the real case) or C" (the complex 
case) with scalar product (, ). Let F = {fir ' ' i fm} be a spanning set of m vectors in 
H . In finite dimension (as it is the case here) such a set forms a frame. In the infinite 
dimensional case, the concept of frame involves a stronger property than completness (see 
for instance pJ3|). We review additional terminology and properties which remain still true 
in the infinite dimensional setting. The set F is frame if and only if there are two positive 
constants Q < A < B < oo (called frame bounds) so that 

m 

AM^ <Y,\{^Jk)? < B\\xf 

k=l 

When we can choose A = B the frame is said tight. For A = B = 1 the frame is called 
Parseval. A set of vectors F of the n-dimensional Hilbert space H is said to have full spark 
if any subset of n vectors is linearly independent. 

For a vector x & H, the collection of coefficients {{x,fk) , 1 < A; < m} represents the 
analysis of vector x given by the frame F. In H we consider the following equivalence relation: 

(2.1) x,y & H , X ^ y iS y = zx for some scalar z with |z| = 1 

Let H = H / ~ be the set of classes of equivalence induced by this relation. Thus x = {x, —x} 
in the real case (when H = M"), and x = {e*°a;, < a < 27r} in the complex case (when 
H = C^). The analysis map induces the following nonlinear map 

(2.2) ip:H^{R+r, ifix) = {\{xjk)\'')i<k<m 

where = {x , x G M , x > 0} is the set of nonnegative real numbers. In previous papers 
[1] we studied when the nonlinear map (p is injective. We review these results below. In this 
paper we describe an algorithm to solve the equation 

(2.3) y = ifix) 

and then we study its performance in the presence of additive white Gaussian noise when 
the model becomes 

(2.4) y = if{x) + z/ , z/ ~ N(0, a^) 

We shall derive the Cramer-Rao Lower Bound (CRLB) for this model and compare its 
performance to this bound. 
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2.1. Existing Results. We revise now existing results on injectivity of the nonlinear map 
Lf. In general a subset Z of a topological space is said generic if its open interior is dense. 
However in the following statements, the term generic refers to Zarisky topology: a set 
Z C K"^*" = K" X ■ ■ ■ X is said generic if Z is dense in K"'^'" and its complement is a 
finite union of zero sets of polynomials in nm variables with coefficients in the field K (here 
K = Mor K = C). 

Theorem 2.1 ([4jTh.2.8). In the real case when H = the following are equivalent: 

(1) The nonlinear map ip is injective; 

(2) For any disjoint partition of the frame set ¥ = Fi UF2, either Fi spans H or ¥2 spans 
H. 

Corollary 2.2 (l20]Th.I;[l]Th.2.2,Prop.2.5,Cor.2.6). The following hold true m the real case 
H = M'^.- 

(1) If if is injective then m > 2n — 1; 

(2) If m <2n — 2 then ip cannot be injective; 

(3) If m = 2n — 1 then (p is injective if an only if ¥ is full spark; 

(4) If m > 2n — 1 and ¥ is full spark then the map if is injective; 

(5) If m > 2n — 1 then for a generic frame ¥ the map (p is injective. 

Theorem 2.3 ([20jTh.II;[l]Th.3.3). In the complex case when if = C" the following hold 
true: 

(1) If m > 4n — 2 then for a generic frame ¥ the map (p is injective; 

(2) // (f is injective then m > 3n — 2; 

(3) If m < 3n — 3 then the map ip cannot be injective. 



2.2. New Injectivity Results. We obtain equivalent conditions to (2) in Theorem 2.1 In 
the real case these new conditions are equivalent to (p being injective; in the complex case 
they are only necessary condition for injectivity. 

Theorem 2.4. Given a m-set of vectors ¥ = {/i, ■ ■ ■ , /^j C H the following conditions are 
equivalent: 

(1) For any disjoint partition of the frame set ¥ = Fi UF2, either ¥1 spans H or ¥2 spans 
H; 

(2) For any two vectors x,y ^ H if n ^ and y ^ then 



fc=i 

(3) There is a positive real constant oq > so that for all x,y E H, 

m 

(2.5) 5^|(x,/,)ri(l/,/.)r>ao||a;f||yf 



k=l 
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(4) There is a positive real constant ao > so that for all x E H, 

m 

(2.6) R{x) ■.= Y,\{xJk)\\-Jk)fk > aol 

k=l 

where the inequality is in the sense of quadratic forms. 
Remark 2.5. The constants in (3) and (4) above are the same (hence the same notation). 
Proof 

(1) =^ (2). We prove by contradiction: -'(2) =^ ""(l)- Assume there are x,y G H, 
x ^ 0, y ^ so that j:T=i\{^Jk)\^\{yJk)\^ = 0- Then {xjk){y,fk) = for all k. Let 

= {k , {x,fk) = 0}, and set Fi = {/^ , k e 4}. Let ly = {!,■■■ ,m}\ 4, and set 
F2 = F \ Fi. Since x is orthogonal to Fi it follows that Fi cannot span the whole H; 
similarly F2 cannot span H because y is orthogonal to all F2. This contradicts (1). 

(2) =^ (3). The unit sphere 5*1 (if) is compact in H and so is Si{H) x Si{H). Since the 
map 

(x,i/)^^i(x,/,)n(y,/,)r 

k=l 

is continuous, it follows 

m 

(2.7) ao := min^^^y)^s^{H)xs^iH)^\{x, fk)\'^\{y, fk)\^ > 0. 

k=l 

By homogeneity for any x,yEH,x^O,yj^Owe obtain: 

m m 
k=l k=l " " 



hence (2.6). If either x = or 7/ = then (2.6) holds true. 

(3) =^ (4). Follows immediately by definition of quadratic forms. 

(4) =^ (1). We prove by contradiction: -i(l) =^ "1(4). If there is a partition F = Fi U F2 
so that neither Fi spans H nor F2 spans H, then there are two non-zero vectors x,y & H so 
that a; _L Fi and y _L F2. Thus (x, fk){y, fk) = for all k. In turn this means y G ker{R{x)) 
which contradicts (2.6). □. 

Note the proof of this result produced the following condition equivalent to negating any 
of the statements of Theorem |2.4 There are two non-zero vectors x,y & H and a subset 
F' C F so that {x, /) = for all / G F', and {y, /) = for all / G F \ F'. Then one can 
immediatly check that x + y and x — y are two non-equivalent vectors in H with respect to 
the (real or complex) equivalence relation ~, and yet (p{x + y) = <f{x — y); hence cannot 
be injective. We thus obtained 

Corollary 2.6. (1) When H = M", (f is injective if and only if any (and hence all) of 
the conditions of Theorem \2.4\ is satisfied; 



(2) When H = , if (p is injective then all conditions of Theorem 2.4 must hold. 
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3. Our approach: Regularized Iterative Least-Square Optimization 



Consider the additive noise model in (2.4). Our data is the vector y G M™". Our goal is to 
find a X G -ff that minimizes \\y — V'(a;)||, where we use the Euclidian norm. As discussed also 
in section |4| the least-square error minimizer represents the Maximum Likelihood Estimator 
(MLE) when the noise is Gaussian. In this section we discuss an optimization algorithm for 
this criterion. Consider the following function 

(3.8) J : if X X R+ X R+ -> M+ 

m 

IVk - {u, fk){fk,v)\ + \\\u\\ +fi\\u-v\\ +A||t;|| . 

k=l 

Our goal is to minimize \\y — f{u)\f = J{u,u,0, fi) over u, for some (and hence any) value 
/i G M"*". Our strategy is the following iterative process: 

(3.9) x*^"*^ = argminuJ{u, , \t, Ht) 
for some initialization and policy {Xt)t>o a^^d {^t)t>o- 

3.1. Initialization. Consider the regularized least-square problem: 

miriuJ^u, u, A, 0) = minu\\y ~ (p{u)\f + 2X\\u\f 
Note the following relation 

m m 

J(«,«,A,0) = Ibf + 2A||Mf -2^|/,|(M,/,)|2 + ^|(M,/fc)|^ 

k=l k=l 
m 

(3.10) = \\yf + 2{{XI-Q)u,u) + Y,\{u,h)\' 



where 

m 

(3.11) Q = 5Zyfc(-,/A.)/fc- 



k=l 



For A > HQ II the optimal solution is -u = 0. Note that if Q < as a quadratic form then 
the optimal solution of mm„||?/ — </)(«) ||^ is m = 0. Consequently we assume the largest 
eigenvalue of Q is positive. As A decreases the optimizer remains small. Hence we can 
neglect the forth order term in u in the expansion above and obtain: 

J{u, u, A, 0) ^ \\yf + 2((A/ - Q)u, u) 

Thus the critical value of A for which we may get a nonzero solution is A = maxeig{Q) 
is the maximum eigenvalue of Q. Let us denote by ei this (positive) eigenvalue and Vi its 
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associated normalized eigenvector. This sugge sts to initialize A = aei for some a < 1 and 
= /3vi, for some scalar (3. Substituting into (3.10) we obtain 

J{/3v,, f3v,, «ei, 0) = \\yf - 2(1 - a)e,P' + (J^ | {v^, h) 

k=l 

For fixed a, the minimum over /3 is achieved at 



(3.12) /3o 



a)ei 



The parameter /i controls the step size at each iteration. The larger the value the smaller the 
step. On the other hand, a small value of this parameter may produce an unstable behavior 
of the iterates. In our implementation we use the same initial value for both A and /x: 

(3.13) /io = Ao = cxei 



3.2. Iterations. Optimization problem (3.9) admits a closed form solution. Specifically, 
expanding the quadratic in u we obtain 

J{u,x\Xt,fj.t) = {{Rt + Xt + fJ't)u,u) - {u,{Q + fit)x^) - {{Q + fJ't)x\u) + 
+ \\yf + {Xt + i^t)\\x'\f 

where 

m 

(3.14) Rt = J2\('''^f^)\'(-^f'^)f^ 

k=l 



and Q is defined in (3.11). We obtain that x*"*"^ satisfies the following linear equation 
(3.15) {Rt + \t + lXt)x'^' = iQ + l^t)x' 

Note the quadratic form on the left hand side is bounded below by 

Rt + Xt + fJ't> ao||a;*||^ + Xt + fit > 



where Oq is given by (2.7). 

3.3. Convergence. Denote jt = J{x^'^^,x^, Xt, fit) = minuJ^u, x^ , Xt, fit) ■ We have the fol- 
lowing general result: 

Theorem 3.1. Assume Xq > Xi > ■ ■ ■ > Xt > ■ ■ ■ and fiQ > fii > ■■■ > fit > ■■■■ 

Then for any initialization x° the sequence {jt)t>o is monotonically decreasing and therefore 
convergent. 

This theorem follows immediately from the following lemma: 
Lemma 3.2. Assume Xt > At+i and fit > fit+i, Then jt > jt+i- 
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Proof First remark the symmetry 



(3.16) 

Then we have: 



J{u, V, A, fi) = J{v, u, A, /i). 



J(a;*^^, x*^"*^, Af+i, < J{x^,x^^^ , Xf+i, fit+i] 



jt+i = J[x ' ,x 

< J{x\ Ai, /ii) = J(x*+\ x\ At, /it) = jt 

This concludes the proof of the lemma. □ 

3.4. First Algorithm. We are now ready to state the first optimization algorithm: 

Input data: Frame set F = {fk , 1 < < m}, measurements y = {yk)t=i^ initialization 

parameter a, stopping criterion e, or maximum number of iterations Tmax- 

Initialization: Compute matrix Q in (3.11 ) and its principal eigenvalue ei and eigenvector 

vi. Compute (3o in (3.12). Set t = and 



X 



Iterate. Repeat: 



(1) Comp ute R t given by (3.14); 

(2) Solve ^A5\) for x'+^; 

(3) Update At, fit using a preset or adaptive policy (more details are provided in section 

(4) Compute jt = J(x*^^, x*. At, fit) and increment t = t + 1; 
Until t > Tmax, or jt-i - jt < e, or At < e. 

Outputs: Estimated signal x*, criterion jt, error \\y — V5(x*)|| . 

3.5. Second Algorithm. Results of numerical simulations (see sectionjs]) suggest the adap- 
tation of A and fi is too agressive. Instead of running the algorithm until A < le — 8 (a small 
value), we implemented a second algorithm where we track the mean-square error: 



L{x') = J2\yk-\{^\h)\ 



2|2 



k=l 



initialization 



and return the minimum value. We thus obtain a second algorithm: 
Input data: Frame set F = {fk , 1 < k < m}, measurements y = {y, 

parameter a, stopping criterion e, or maximum number of iterations Tmax- 

Initialization: Compute matrix Q in (3.11 ) and its principal eigenvalue ei and eigenvector 

vi. Compute (3o in (3.12). Set t = and 

Xo = fio = otei , x° = PoVi , L 



optim 



^ ^ l^fc \{^ , fk)\ I ; ^ optim 



X 



k=l 



Iterate. Repeat: 

(1) Comp ute R t given by (3.14); 

(2) Solve ([3l5|) for x*+i; 
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(3) Update At, [it using a preset or adaptive policy (more details are provided in section 

i; 

(4) Compute jt = J(x*"'"^, x*, Xt, fit)', 

(5) Compute = ELi IVk - \{x'^\ fk)\r; 

(6) If Lt^i < Loptim then LopUm = -^t+i aud XopUm = 2: ; 

(7) increment t = t + 1; 

Until t > T^ax, or - < e, or A* < e. 

Outputs: Estimated signal XopUm, criterion jf, error LopUm = \\y - <^{xoptim)\f ■ 

4. The Cramer-Rao Lower Bounds 



Consider the noisy measurement model (2.4), y = ^{x) + z/, with v ~ N(0,(T ). Fix a 
direction in if, say e. We make the following two assumptions regarding the unknown signal 
x: (1) We assume x is not orthogonal to e, that is (x, e) 7^ 0; (2) We assume we are given 
the sign of this scalar product; say (x, e) > 0. These two assumptions allow us to resolve the 
global sign ambiguity. Thus x E S d H where 5" is a half-space of H. Since it is a convex 
cone we can compute expectations of random variables defined in S. The likelihood function 
is given by 

The Fisher information matrix I(x) is given by 

where the expectation is over the noise process, for fixed x. 

In the following we perform the computations in the real case H = M."'. For ease of 
notation we assume the canonical basis of M" and the lower index (or second index) denotes 
the coordinate with respect to this basis; for instance Xk and fi^k denote the k^^ coordinate 
of X and respectively. 

^^^^ - - - ^) 4 E « (K^. /.>i- ^0 

^ 1=1 1=1 

Now use 'E[yi] = We thus obtain 

(4.17) I(x) = -,f2\{xJi)\'fiff = -^R{x) 

1=1 



where R{x) denotes the quadratic form introduced in (2.6). Now we are ready to state the 
first lower bound result (see e.g. [2Z] Theorem 3.2). 



RECONSTRUCTION OF SIGNALS FROM MAGNITUDES OF REDUNDANT REPRESENTATIONS 9 



Theorem 4.1. In the real case H = M", the Fisher information matrix for model (2.4) is 
given by I{x) in (^-IT)- Consequently the covariance matrix of any uniahsed estimator $(?/) 
for X is hounded below by the Cramer-Rao lower bound as follows 

2 

(4.18) Var[^{y)] > CRLB{x) := (lix))-^ = ^(i?(x))-^ 

Furthermore the conditional mean-square error of any unbiased estimator $(?/) is given by 

2 

(4.19) E[||<l>(?/) - xf\x] > ^trace{{R{x))-^}. 

When signal x is random and drawn from x ~ N(0, /), the mean-square error of the unbiased 
estimator $(y) is bounded below by 

2 

(4.20) E[||$(?/) - xf] > ^trace{E[(i?(x))-i]}. 



Remark 4.2. Corollary \27b\ implies that when ip is injective the Fisher information matrix 
is invertible, hence a bounded CRLB, and the signal is identifiable in S (up to a global phase 
factor). 

We derive now a different lower bound for a modified estimation problem. Let us denote 
by X = xx'^ and Fk = fkfk the rank-1 operators associated to vectors x and fk respectively. 
Note \ {x, fk)\'^ = trace{XFk}. Hence 

yk = trace{XFk} + Vk ■, 1 < k < m 

We would like to obtain a lower bound on conditional mean-square error of an unbiased 
estimator of the rank-1 matrix X. A naive computation of the Fisher information associated 
to X in the linear model above would produce a singular matrix whenever m < "-^"-^^^ (the 
reason being the fact that a general symmetric X is not identifiable merely from m < ^^'^^^^ 
measurements). Instead the bound should be derived under the additional hypothesis that 

2 

X has rank one. We obtain such a bound using a modified CRLB. Let (7 : S* — ?► M" be the 
vector valued map 

gix)= ( x,x, )i<,<^.<„ 

of components. Let "^{y) denote any unbiased estimator of the rank-1 matrix X. Then 
(see equation (3.30) in [27j ) 



(4.21) Co„(*fe))>^g|jI 
Taking trace on both sides we get 



99 \ ^-1 f dg 



T 



Emiy)-xx^\f\x] > trace{r' (^^J } 
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Let H = Then we have 



dx ' 



z,j=l 



Thus we obtained 



Theorem 4.3. The conditional mean-square error of any unbiased estimator of the rank-1 
matrix X = xx'^ is bounded below by 

2 

(4.22) E[||^(?/) - xx^f |x] > y {\\xftrace{R-^} + x^R'^x) . 

Consider now the case of the Maximum Likehhood Estimator (MLE) whose optimization 
problem was considered in the previous section. For model (2.4) this takes the form of 

(4.23) ^MLE{y) = argmiUuWfiu) - y\f = argminuJ{u, u; 0, 0) 

The MLE computes the global minimum in the optimization problem above. Assume that 
^MLE selects the closest global minimum to x. We want to estimate lower bounds on the 
MLE performance so that we can benchmark performance of any optimization algorithm 
against these bounds. 

First we estimate the bias of the MLE estimator in the asymptotic limit a — )■ 0. The 
estimator must satify the MLE equation 

V(/oc/L(x))U.=$,,^^(j^) = 

which turns into 

m 

J2 i^MLEiy), fk) ili^MLEiy), fk)\^ - yk) fk = 

k=l 

Denote v = ^MLE^y) — x. The bias is given by Bias{x) = K[v\x]. Assymptotically we 
can assume ||f || is small with high probability. We shall expand the products in the above 
equation taking into account only the first terms in v. 

m 

J2({x, fk) + {v, fk)mx, fk) + {v, fk)\' - yk)fk = 

k=l 

Expanding the products and neglecting higher order terms in v we obtain: 

(rn \ m 

$^(3|(a;, fk)\' - yk)fkfl V + J2{\{x, - y,){x, f,)f, = 
k=l / k=l 
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Note Uk = \ {x, + z/fc. Let us denote 

m m 

N = Y.UkfkfI , R = R{x) = Y,\{x,fk)\'fkfI 

k=l k=l 

The equation that v satisfies becomes (2i? — N)v — Nx = 0. Therefore 
V = {2R- N)-^Nx Bias{x) =E[{2R~ N)-^N]x. 



For fixed x 0, due to the lower bound in (2.6) we obtain with high probabihty ||A^|| < 
2eigmin{R), where eig min{R) denotes the smallest eigenvalue of R. Note eigmin{R) = 
> ctolFlr by (2.7). Then using Neumann's series expansion we get 

^(3| (x, - yk)fkfi] = {2R - N)-' = ^R-' + -R-'NR-^ + 0{N') 



^k=l 

Thus we obtain 



v = {2R- N)~^Nx = ^R-^Nx + ^R-^NR'^Nx + 0{N^) 

Note also the similarity criterion in expansion above is related to . which is of the 

order (t||-R~^||. Since all odd moments of Gaussian random variables vanish we obtain 
E[(i?-iiV)2p+i] = 0. Hence 

Bias{x) = E[v] = ^E[R'^NR~^Nx] +E[0{N^)] = 

2 ™ 

(4.25) = ^J2{x,fk){R''fk,fk)R-'fk + 0{ia\\R-'\\y) 

k=l 

The leading term in bias has the form 

(4.26) 5ms°(x) = ^6 , where 6 = {^Jk){{Rix))-' fk, fk){Rix)y' fk 

k=l 

Note the dependence on x is highly nonlinear. We would like next to obtain the modified 
CRLB for MLE taking into account its bias. We need to estimate the first derivatives of 



Bias{x) with respect to x, ^liEl^fiElh. _ Again we shall derive the asymptotic approximation 
of this matrix: 

_ 4 d{Bias%x))j _ d6j{x) 
cr^ dxi dxi 

The key relation to use is 

— {R{x))-^h = -R-'g^R-'fk = -2j2{x,U)U,i{R-'fk,U)R-'U 
' ' p=i 



which comes from R{x){R{x)) ^fk = fk by differentiating with respect to xi, and from (2.6). 
After some straightforward but tedious algebra we obtain 
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(4.27) A = 

k=l 
m 

- 2 {x,h){x,fp){R-'fk,h){R-'fk,fp)R-'fpff 

k,p=l 



- 2 J2 {xjk){x,mR-'fpjk)\'R-'fkfi 



ip 

k,p=l 



Now we can compute the modified Cramer-Rao lower bound for tlie MLE estimator (see 
e.g. [27] Equation (3.30)). 



Theorem 4.4. The MLE estimator (4-23) is biased. Its expectation admits the following 
asymptotic approximation 

2 

(4.28) E[^MLE{y)] =x + ^5 + Oi{cj\\R-'\\Y). 

Its covariance matrix is bounded below by 
(4.29) 

Cov[<i>MLE{y)] > = ^R-'+^QiR-'A^+AR-')+omR-^r) 

where I is the identity matrix. Furthermore, the conditional mean-square error is bounded 
below by 



(4.30) 



E[\\^MLE{y) - xf\x] = \\Bias{x)\f + trace{Cov[^MLE{y)]} 



> —trace{R-^} + — (||5|r + 2trace{AR-H) + 0((a\\R-^\ 
4 16 



Here we used the notation R = R{x) = EfcLi /fc)! /fc/fc ; and 6, A given by (4.26), (4-21}. 



5. Numerical Analysis 

In this section we present numerical simulations for the algorithms presented in this paper. 

We generated random frames or redundancy 3, that is m = 3n, as well as random signals 
X. All these vectors (frame and signal) are drawn from N(0,/). We set the first component 
of X positive, and so we decided the global sign after reconstruction. To the magnitude 
square of signal coefficients v?(x) we added Gaussian noise of variance cr^ to achieve a fixed 
signal-to-noise-ratio defined as 

\lr fi.)\^ 

SNR = ^^=^ r^^^ , SNRdB = 10 login (^iVi?) [dB] 
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Note the similarity criterion used in asymptotic expansions (4.28) and (4.29) is of the same 
order as a 11 i?^ 11 ~ ^gjyj^ (^P multiphcative constants). We used 11 values of SNRdB in 
lOdB increments from -20dB to +80dB. 

For the first algorithm, results are averaged over 100 noise realizations. In each instance 
of the algorithm we initialized Aq, /xo as described in subsection 3^ with a = 0.9. At each 
iteration Xt+i = Aj/1.05. We run the algorithm for at least 100 steps, or until At gets below 
10~®. The parameter/it is adapted as follows: = Tnax{l, Aj). 

We include results for n = 10, n = 50 and n = 100. Figure [T] includes the conditional mean- 
square error averaged over 100 noise realizations, and the lower bounds: the unbiased CRLB 
(4.18) and the MLE adapted CRLB (4.30). Note the two lower bounds are indistingueshable 
for SNR > 20dB. For low SNR, when the two bounds differ significantly, the approximation 
(4.24) is no longer valid. Hence the bound would be different as well. In general we cannot 
obtaion a closed form solution for the new bound. 

In Figure |2] we plot the bias and variance components of the mean-square error for the 
same results in Figure [T] Note the bias is relatively small. The bulk of mean-square error is 
due to estimation variance. 

Figure |3] contains the average number of iterations for each of these cases. The algorithm 
runs for about 530-660 steps. 

For the second algorithm we repeated the same cases [n = 10, 50, 100) and same levels of 
SNR, but we average over 1000 noise realizations. We present the mean-square error in two 
cases: in Figure |4] the case of fixed sign as discussed earlier (first component of x is positive); 
in Figure [5] the case of a sign oracle, when the global sign is chosen as given by the minimum 

?77.Z7z( II X 3^ optim II 5 ll"'' optim || ) • 



6. Conclusions 

Novel necessary conditions for signal reconstruction from magnitudes of frame coefficients 
have been presented. These conditions are also sufficient in the real case. The least-square 
criterion has been analyzed, and two algorithms have been proposed to optimize this cri- 
terion. Performance of the second algorithm presented in this paper is remarkably close to 
the theoretical lower bound given by the Cramer-Rao inequality. In fact for low SNR its 
performance is better than the asymptotic approximation given by the modified CRLB. 
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Figure 1. Conditional Mean-Square Error and CRLB bounds for n 
(top plot), n = 50 (middle plot), and n = 100 (bottom plot). 
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Figure 2. Bias and Variance components of conditional Mean-Square Error 
and CRLB bounds for n = 10 (top plot), n = 50 (middle plot), and n = 100 
(bottom plot). 
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Figure 3. Average number of iterations vs. SNR for n = 10 (top plot), 
n = 50 (middle plot), and n = 100 (bottom plot). 
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Figure 4. Conditional Mean-Square Error and CRLB bounds for n = 10 
(top plot), n = 50 (middle plot), and n = 100 (bottom plot) for fixed sign. 
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Mean square error of the sign oracle estimator 
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Figure 5. Conditional Mean-Square Error and CRLB bounds for n = 
(top plot), n = 50 (middle plot), and n = 100 (bottom plot) for global si 
oracle. 



